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Abstract. We present techniques to perturb, measure and model the ion velocity distribution in an 
ultracold neutral plasma produced by photoionization of strontium atoms. By optical pumping with 
circularly polarized light we promote ions with certain velocities to a different spin ground state, and 
probe the resulting perturbed velocity distribution through laser-induced fluorescence spectroscopy. 
We discuss various approaches to extract the velocity distribution from our measured spectra, and 
assess their quality through comparisons with molecular dynamic simulations. 
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INTRODUCTION 

The ability to modify and probe the velocity distribution of particles in plasmas is 
valuable for studying collective modes [1], transport [2], and thermalization rates in 
plasmas [3] . Here, we describe the combined application of optical pumping and laser- 
induced fluorescence (LIF) Doppler spectroscopy to create and characterize perturbed 
ion velocity distributions in ultracold plasmas. Optical pumping [4] between different 
hyperfine states provides a powerful method for tagging certain particles in equilibrium 
systems. Exploiting the Doppler effect, the spin-state modification can be done in a 
velocity-dependent manner to modify the velocity distribution of a given spin state. 
On the other hand, LIF spectroscopy is a well established tool for measuring velocity 
distributions [5] and has been used to probe pure ion plasmas [1] or plasmas created with 
short pulse lasers [6]. 

Ultracold plasmas (UCPs), produced by photo-ionization of laser-cooled atoms [7] 
or cold molecular beams [8] provide a well controlled laboratory to study various 
plasma physics phenomena, such as collective waves [9, 10, 11, 12, 13, 14, 15], plasma 
expansion into vacuum [16, 17, 18, 19, 20, 21, 22, 23], correlation effects [24, 25, 26, 
27, 28, 29, 30, 31, 32], recombination to form neutral atoms [33, 34, 35, 36, 37, 38], and 
plasma instabilities [39, 40]. Due to their very low temperatures UCPs also realize an 
interesting parameter regime [41, 42, 43] in which the ionic plasma component can 
be strongly coupled. Consequently, the present work opens up new studies of non- 
equilibrium plasma dynamics in the strong coupling regime. 
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FIGURE 1. Schematics of the experimental setup. Neutral strontium atoms from a heated reservoir are 
Zeeman-slowed before entering the trapping region. The magneto-optical trap consists of a pair of anti- 
Helmholtz magnetic coils and six laser beams. Pi trapped atoms are then ionized by the photoionizing 
laser. The fluorescence probe beam propagates in a direction that is perpendicular to the imaging axis 
and CCD camera. The complementary absorption probe beam passes through the plasma and falls on the 
camera. Adapted from [28]. 



CREATION OF ULTRACOLD NEUTRAL PLASMAS 



We create an ultracold neutral plasma through photoionization of laser-cooled strontium 
atoms in a magneto-optical trap (MOT) [28, 44, 45]. Figure 1 shows our experimental 
setup. The MOT operates on the dipole- allowed So — 1 Pi transition of 88 Sr at 461 
nm, with transition linewidth, y/(2n) = 30.5 MHz [46]. The laser-cooled atom cloud 
is characterized by a temperature of ~ 10 mK and has a spherically symmetric Gaussian 
density distribution, n(r) = rcoexp(— r 2 /2o 2 ), with o ~ 0.6 mm and hq w 6 x 10 10 cm~ 3 . 
The number of trapped atoms is typically 2 x 10 8 . 

The atoms are ionized via two-photon ionization by two temporally and spatially 
overlapping, retro-reflected ~10 ns laser pulses: the first is obtained from a pulse- 
amplified laser beam tuned to the cooling transition of the atoms ( So — 1 Pi) at 461 
nm and the second one derives from a pulsed dye laser tuned just above the ionization 
continuum at ~ 412 nm (Figure 2 A). 

In this way we ionize ~30-70% of the atoms. The plasma inherits its density distri- 
bution from the neutral atoms, resulting in peak electron and ion densities as high as 
noe ~ n 0i ~ 4.2 x 10 10 cm~ 3 . The remaining ground state atoms have no effect on the 
subsequent plasma dynamics, due to the short time scale of the experiment and the small 
neutral-ion collision cross- sections. 

As a result of the much lighter mass of the electrons (compared to the heavy Sr ions), 
most of the excess energy from the photoionizing beam is acquired by the electrons, 
while the ions' kinetic energy remains similar to that of the neutral atoms in the MOT 
[41]. By tuning the wavelength of the pulsed dye laser, we vary the initial electron kinetic 
energies typically between 1 and 1000 K. The ion kinetic energy remains on the order of 
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FIGURE 2. Relevant atomic and ionic energy levels of strontium, with corresponding decay rates. (A) 
Neutral atoms are laser-cooled and trapped in a MOT operating on the 1 So — X Pi transition at 460.9 nm 
[45]. Atoms are excited to the l P\ level by a pulse-amplified laser and ionized by pulsed dye laser at 



~ 412 nm. (B) Imaging of the ions is done on the 2 S 1/2 — 2 Pi/2 transition at 421.7 nm. 2 P\/2 ions decay to 
the 2 D3/2 state 7% of the time, after which they cease to interact with the probe beam. The intensity and 
duration of the 422 nm light is sufficiently low to avoid optical pumping to the metastable 2 Dj,/2 state. 



millikelvin, close to the kinetic energy of neutral atoms in the MOT. However, disorder- 
induced heating [24, 29, 30] raises the temperature of the cold ions to approximately 1 
K on the timescale of the inverse ion plasma oscillation frequency (~ 500 ns). 



OPTICAL MODIFICATION OF VELOCITY DISTRIBUTIONS 

The internal structure of the strontium ions relevant to our experiment can be represented 
by the simplified four level scheme depicted in Figure 3. Both the ground state (5 2 S 1 / 2 ) 
and the excited state (5 2 Pi/ 2 ) are doubly degenerate, where the total angular momentum 
projection can take on the values mj = +1/2 and mj = —1/2. The two ground states 
of the 5 2 S 1 / 2 level will be denoted by 1 (mj = —1/2) and 3 (my = +1/2), and the two 
excited states of the 5 ¥1/2 level will be denoted by 2 (mj = +1/2) and 4 (mj = —1/2). 

The transitions between level 1 and 2 and between level 3 and 4 are driven by lasers 
with circularly polarized light, 0+ and o respectively, with frequencies C0y and Rabi 
frequencies 

where Eq is the field amplitude of the laser. The corresponding frequency detuning for 
each transition is Ay = C0y — co, where ha = AE is the energy difference between ground 
and excited state. 
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FIGURE 3. Simplified energy level schematic for the Sr + , where dashed lines indicate decay paths 
while solid lines indicate excitation paths. The mj = — 1/2 (mj = +1 /2) spin ground level is coupled to 
the higher energy mj = +1 /2 (mj = —1/2) via a + (ex - ) circularly polarized laser beam, from where it 
can decay back to either the mj = — 1/2 or mj = +1/2. n ; - labels the population of a particular spin state. 



Ions in the excited states (2 and 4) decay to the ground states (1 and 3) with decay 
rates 



where is the dipol matrix element of the optical transition between level i and j, X 
is the wavelength of the transition between ground and excited state, £q is the electric 
constant and h is the reduced Planck constant. For the total decay rate of levels 2 and 4 
we have T = T\2 + ^23 — + ^4, respectively. 
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Employing the rotating wave approximation one obtains the following optical Bloch 
equations for the density matrix of the four-level system 

Pi l = i -y (P21 - P12) + Ti2 P22 + Ti4 p44 (3a) 

P22 = ~i -y (P21 - P12) - r P22 (3b) 

P33 = i -y- (P43 - P34) + T 2 3 P22 + T 34 p44 (3c) 

p44 = -i y^ (P43 - P34) - r p44 (3d) 

Q12 r 

P12 = -i -y- (Ph _ P22) - ( 2 + i A12) P12 (3e) 

P13 = i P23 - 1 P14 (3f) 

r 

P14 = i P24 - i -p- P13 - ( 2 + i A34) Pl4 (3g) 
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P23 = i Pl3 - i -7T P24 - ( 2 - i A12) p23 (3h) 

P24 = i -y Pl4 - i -y P23 - (r - i (A i2 - A34) )p 24 (3i) 

q 34 r 

P34 = ~i -y (P33 - P44) - ( 2 + i A34) P34- (3j) 

The diagonal elements p ;; are the populations of the four levels i = 1,2,3,4 and pij, i 7^ j 
represent the coherences of the corresponding transitions. 

At long times t > T 1 one can adiabatically eliminate the coherences, i.e. set p,j = 
for i 7^ j, to rewrite Equations (3) in terms of a simple set of rate equations 

/i(v) = +Rn(v)[f 2 (v)-Mv)} +r 12 / 2 (v) + r 14 / 4 (v), (4a) 

/ 2 (v) = -/?i 2 (v)[/ 2 (v) -/i(v)] - (r 12 +r 23 )/ 2 (v), (4b) 

/ 3 (v) = +/?34(V)[/ 4 (V) -/ 3 (V)] +r 23 / 2 (v) +r 34 / 4 (v), (4C) 

/ 4 (v) = -i? 34 (v)[/ 4 (v) -/ 3 (v)] - (r 14 +r 34 )/ 4 (v), (4d) 

for the velocity distributions /, = pa of ions in a given spin state i. The resulting pumping 
rates have a simple Lorentzian shape, 

*^=iT[2OT (5) 

As in the experiment two counterpropagating laser beams with same detuning A pump and 
wavenumber k = 2n/X are used, the Doppler shift of an ion moving with velocity v in 
the directions of the lasers beams is accounted for by replacing Ai 2 = A pump + kv and 
A 34 = A pump — kv in Equations (3) and (5). 

The simplified treatment of the rate equations (4) shows that the laser pumping 
transfers spin population between level 1 and 3 in the domain around the resonant 



N(t) 





A 

fm 


i 


* 
* 
■ 




/ * 


m 
m 
m 

« 






. v(t) 









-v Res + v Res 

FIGURE 4. Cartoon of the expected velocity distribution for the 5 2 Si/2 ('«/ = +1/2), «3, with both 
counterpropagating pump lasers turned on (d - and The solid line indicates the unperturbed distribu- 
tion, while the dashed line indicate changes induced by optical pumping. A depletion of ion population is 
created as ions that meet the resonance condition VR es = A pump /£ with the <7~ laser while an enhancement 
in ion population appears as ions from (mj — —1/2) 5 2 Si/2 ( n i) level are optically pumped by the a + 
laser. 



velocities v res = ±A pump /fc, creating thereby an asymmetric velocity distribution as 
sketched in Figure 4. 

Our experimental scheme to perform the described optical pumping and the resulting 
velocity distributions is shown in Figure 5. The two circularly polarized pump beams, 
with the same frequency yet opposite polarization and direction, interact with the plasma 
for a time of ~ 100 ns. For the strontium ion, the S 1/2 ~~ P1/2 transition has a wavelength 
A = 421.7 nm, and a total decay rate Y / (2n) = 21 MHz. The a + beam has an intensity 
of = 334mW/cm 2 with beam waist of 1.18 mm, resulting in Q12/ (2n) = 36.6 MHz. 
On the other hand, the pumping beam has an intensity of I a _ = 315 mW/cm 2 with 
beam waist of 1 . 12 mm, resulting in ^34/ (2n) = 35.0 MHz. Furthermore, A pump / (2n) = 
—20 MHz. For these parameters the resonant velocity |v res | = 8.4 m/s is well within 
the thermal velocity ^/k^T jm = 15 m/s allowing for an efficient modification of the 
velocity distribution. 

MEASUREMENT OF PERTURBED VELOCITY SPECTRA 

At an adjustable time after optical pumping, a third, less intense o circularly polarized 
probe beam, near resonant to the transition between level 3 and 4, propagates through 
the plasma (see Figure. 5). This probe beam has an intensity I pro t, e = 73 mW/cm 2 with 
beam waists of w z = 0.625 mm and w x = 4.68 mm. Laser-induced fluorescence in a 
perpendicular direction is then captured by a camera. The signal is spatially-resolved, 
so we can analyze the signal from different regions near the center of the plasma that 
have little variation in density and average plasma expansion velocity. By scanning the 
frequency detuning A of the probe beam, we obtain a spectrum 5(A) that contains 
information on the ion velocity distribution. The measured fluorescence spectrum is 
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FIGURE 5. Experimental setup for optical modification of the ion velocity distribution in an UCP. The 
plasma cloud is illuminated by both the pump and probe beams originating from a 422 nm laser, that is 
near resonant with the S1/2 — P1/2 transition of Sr + . The pump beams have the same laser frequency 
yet opposite circular polarization and travel in opposite directions. The probe beam has a circular 
polarization and can be scanned to obtain a fluorescence spectrum by capturing the probe-beam induced 
fluorescence light with a CCD-intensified camera perpendicular to the beam direction. 



proportional to the velocity-dependent population of the ground state f$ convolved with 
the Lorentzian excitation profile of the probe beam, i.e. 

5(A) - J f 3 ( v )L{v,A)dv, (6) 

L(V ' A) ~ l+^probe + [2(A-fcv)/r tot ]2' (7) 

where, r tot = r + rj nst is the total linewidth, rj nst / (2n) = 7 MHz is the instrumental laser 
linewidth and so.probe = Ipmbe/hat- The saturation intensity for the Sm — P1/2 transition 
with circular polarized light is / sa t = 57 mW/cm 2 . 

Figure 6 shows ion spectra measured after 100 ns of optical pumping and, for compar- 
ison, ion spectra without optical pumping, for a plasma with ion temperature 7} = 2.3 K, 
density hq = 3. 1 x 10 15 m~ 3 and coupling parameter T, = 1.7. The spectra are scaled by 
the maximum of each curve. We record separately the spectra for three regions around 
plasma center (labeled Region -1,0, 1). Given the higher mean velocity in the outer re- 
gions due to plasma expansion [47], their corresponding spectra are shifted from zero 
(Regions ±1). Notice the depletion in ion population for A/ (2n) ~ —20 MHz and en- 
hancement in ion population for A/ (2n) m 20 MHz in region 0, as discussed above. The 
features here are less sharp, compared to the schematic distribution of Figure 4 due to 
power broadening from the pump beams [44] and the convolution of the velocity dis- 
tribution and Lorentzian excitation profile (Equation 6). As the probe and pump pulses 
originate from the same laser source, the probe pulse has a limited frequency scanning 
range, resulting in the narrow range from about -50 to 50 MHz (see Figure 6), which 
complicates the analysis of the experimental data. 
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FIGURE 6. Ion spectra for 3 regions around plasma center with (red dots) and without (blue dots) pump 
beams for T, ■ = 2.3 K, no = 3.1 x 10 15 m~ 3 and T ; - = 1.7. The spectra were measured 0.1 3 /is after optical 
pumping and are normalized by scaling to the maximum of each curve. To achieve low T,, we allow the 
plasma cloud to expand significantly before applying the pump beams. As a consequence regions other 
than plasma center have significant expansion velocity and display significant shift in their spectra. Notice 
the depletion in ion population for A/(2n) w —20 MHz and enhancement for A/(2n) « +20 MHz. 



EXTRACTING THE VELOCITY DISTRIBUTION 

As described in the previous section, the reconstruction of the underlying velocity 
distribution from the recorded ion spectra is complicated due to experimental limitations 
on the accessible frequency range. This leads to a loss of information contained in the 
wings of the spectrum. To determine possible errors of the reconstruction process, we 
performed numerical simulations of the experiment, which provide spectral information 
over the full frequency range, and can thus be used for accurate comparisons of the 
real and extracted velocity distribution, as obtained from a finite frequency range of the 
simulated spectra. 

The ion dynamics including the collisional redistribution of velocities can be accu- 
rately captured by classical molecular dynamics (MD) simulations. The ions are rep- 
resented by an one-component plasma of N particles in a cubic simulation cell with 
periodic boundary conditions. The mutual interactions of the ions are effectively treated 
by the fast multipole method [48, 49, 50], which permits force calculations for large 
particle numbers with a numerical effort that scales only linearly with the number of 
ions. 

The optical pumping process is described by solving the optical Bloch equations (3) 
along each ion trajectory. The internal state dynamics is included by assigning the four- 
level density matrix pW to each of the q = l,...,N ions. Initially the two ground states 

are equally populated, i. e. p[f = p^f = 0.5, and pj^ = for all remaining matrix el- 
ements. For the duration of the pump pulse we solve the optical Bloch equations (3) 
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FIGURE 7. Comparison of the true average velocity (v) (blue line) to the average velocity as obtained 
from different extraction procedures using the frequency range (— A„,,A m ) of the simulated spectrum. 
Direct integration of the spectrum using Equation (9) and the shifted Maxwellian center velocity v c show 
significant deviations within the experimentally available frequency range A„,/ (2n) sa 50 MHz (shaded 
grey), whereas excellent agreement for the average velocity is obtained from the fit formula (13). 



numerically for each ion individually, taking into account the time-dependent laser de- 
tuning due to the changing ion velocity obtained from the MD simulation. From the 
simulated velocity distributions fi(v,t), we obtain the fluorescence spectra by convolu- 
tion with the Lorentzian emission profile (7) according to Equation (6). 

With the simulated spectra at hand, we can now proceed to analyze possible schemes 
to extract the velocity distribution from the finite frequency range of the experimental 
fluorescence spectra. To assess the quality of the different methods, we use the average 
velocity (v) as a figure of merit. 

Given the entire excitation spectrum, 5(A), the average velocity is related to the 
average frequency shift through the wavenumber of the probe beam, k, according to 

/ AS(A)dA 

(A) = ==5 = k(v) . (8) 

/ S(A)dA 
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FIGURE 8. Fit (green line) of a convoluted Maxwellian (left) and of equation (13) (right) to the 
simulated spectrum (symbols) corresponding to T = 2.8 K, p = 13 x 10 9 cm~ 3 , r = 1.4. In each case 
the fit is only applied in the interval (—50,50) MHz (shaded gray area). 

From a limited portion of the spectrum, — A m < A < A m , if one makes a Taylor expansion 
of the integrand in Equation (6) around v = 0, the relation can be written as 

/ A5(A)JA 

(A>==£ ^K(A m )(v> (9) 

I S(A)dA 

The frequency-interval dependent proportionality constant is given by 

= * f 1 - ] (10) 

\ (l + (2 A„ 7 /r) 2 J arctan (2 A m /T) J 

and recovers Equation (8) in the limit A m — > °°. As shown in Figure 7 this procedure 
yields the correct average velocity for large A m , but significantly underestimates the 
correct value of (v) for the experimentally measurable region (A m ~ 50 MHz). 

An alternative method, which works well for our data, is to assume a specific form for 
the velocity distribution and vary the parameters of the distribution to fit the resulting 
excitation spectrum to the experimental measurements. For an appropriate choice of the 
functional form of the velocity distribution, this allows accurate approximation of the 
full velocity distribution based on finite spectral information. 

A simple ansatz is to assume that the velocity distribution after pumping still retains 
its Maxwellian shape but is only shifted by an amount v c from zero. Figure 8a shows 
the fit of a convoluted, shifted Maxwellian to the simulated spectrum using only the 
experimentally available frequency range (shaded gray area). We find good agreement 
between fit and data in the central, observable part of the spectrum, but slight deviations 
at larger frequencies beyond the experimentally available data. These small discrepan- 
cies can cause considerable deviations of the average velocity as shown Figure 7. 

The situation can be improved by using a fit formula that captures the effects of 
optical pumping and the collisional redistribution of velocities more accurately. In order 



to construct such an expression we start from the rate equations (4) and adiabatically 
eliminate the dynamics of the excited state populations / 2 and f\. Collisional effects are 
accounted for within the relaxation time approximation [51] which gives 



h (v) = (v) h (v) - (v) fi(v)-y 
h (v) = (v) A (v) - (v) f 3 ( v) - r 



/i(v)-^/m(v) 
/3(v)-^/m(v) 



(11a) 
(lib) 



where 7 is an effective relaxation rate and /m (v) = /1 (v) + /3 (v) denotes the Maxwellian 
velocity distribution at a given temperature. Moreover, the transition rates R^, and 7?^ 
for a transition from state 1 to state 3 and vice versa are obtained from the rate equations 
(4) according to 

tfi2(v)r 23 

^->(v) = — r^Ti (12a) 

^ l J r + J?i 2 (v)' 

From the steady state (/i(v) = /3(v) = 0) of these equations we obtain the following 
simple fit formula 

/3(v)=/M(v) ^ g^ 1 f 03) 

l+*l(v res -v) 2 + l+fe 2 (Vres+v) 2 + 

For steady state conditions, the parameters a1.a2.b1 and ^2 can be expressed in terms 
of the parameters of Equations (11) and (12). In most cases, however, the system does 
not reach the steady state for our experimental parameters and optical pumping times. 
Nevertheless, one can apply Equation (13) to fit the measured spectra using a\ , 0,2, b\ , Z?2 
as free parameters. As we will show below, this ansatz indeed provides an excellent de- 
scription of the actual velocity distributions. The parameters a 1,^2 quantify the impor- 
tance of optical pumping relative to collisional redistribution and the parameters b\ , bi 
control the widths of the asymmetric features of the velocity distribution. In the limit 
ai,a 2 3> 1, Equation (13) yields the correct steady state for very strong pumping and 
weak collisions. On the other hand, for a\,a,2 = it recovers the equilibrium Maxwell 
distribution, as expected in the absence of optical pumping. 

As shown in Figure 8b, the frequency spectrum obtained from convolving Equation 
(13) according to Equation (7) yields a much improved fit as compared to the shifted 
Maxwellian, shown in Figure 8a. Moreover it accurately describes the wings of the 
spectrum, even though the fitting has been determined from the finite frequency range 
—50 MHz < A < 50 MHz. Consequently, the extracted average velocity, shown in Figure 
7, is in excellent agreement with the real value of (v) even for a very limited range 
of frequencies available for fitting the fluorescence spectrum. Moreover, the procedure 
described above can also be used to construct reliable fitting functions for more complex 
optical pumping schemes that may be applied in future experiments. 



SUMMARY 



In this work we have presented new techniques to perturb, probe, and model ion velocity 
distributions in ultracold neutral plasmas. We anticipate that the described approach will 
serve as a valuable tool for direct experimental measurements of non-equilibrium plasma 
dynamics in the regime of strong coupling. 
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